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ABSTRACT 



The problem of orbital parameter estimation using angles only observations is 
examined. Direction cosine measurements, obtained from satellite passage of an 
earth-based stationary planar radar beam, are assimilated by an extended Kalman 
filter to improve estimates of a classical orbital element set. Several progressively 
comprehensive orbital motion models are considered and compared. The relative 
effectiveness of these models is illustrated by applying them to actual satellite 
data. 
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I. INTRODUCTION 



Earth-orbiting satellites are catalogued in the form of orbital element sets. For 
maximum usefulness of these catalogs to military and civil space systems, element 
sets need to be updated as frequently and as accurately as possible. A nonlinear 
filter was designed to perform these updates in near-real-time using observations 
of satellites piercing a planar radar beam. 

The filter predicts time-varying orbital parameters according to a non-linear 
orbital dynamics model. This model includes, along with two-body Keplerian 
motion, first order oblateness effects, both secular and periodic, simple decreasing 
exponential atmospheric drag effects, and white noise on both the dynamic model 
and the measured observations. This simple model is found sufficient to track 
several selected satellites, with proper tuning of the filter. It should be noted that 
even this simple non-linear model gives rise to comparatively complex filter 
expressions. 

Applying this filter to the orbit improvement problem helps meet the need for 
continual determination of orbit decay, collision avoidance, satellite maneuvers, 
etc. Current systems use a daily batch least squares differential correction method. 
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II. REVIEW OF LITERATURE 



This research involves the integration of a simple, yet adequate, orbital 
dynamics model with a necessarily suboptimal extended Kalman estimator. This 
combination is applied to a geometrically unique detection scheme to provide near 
real-time improvement of orbital elements for selected satellites. The accuracy 
level under consideration falls under the umbrella of the general perturbations 
problem of orbital mechanics, which is covered by numerous textbooks, a 
sampling of which are referred to here [1-6]. In general, solution procedures 
involve a gravitational potential function (earth oblateness, third body, etc.) 
describing its contribution to the perturbation of a classical two-body Keplerian 
dynamic model. The resultant expressions take the form of variations in some set 
of orbital parameters. The two broad divisions of orbital parameter representation 
in the literature are: 1) some set of 6 orbital elements which describe the orbital 
path as a conic section and its orientation in inertial space, and 2) Cartesian 
position and velocity of the orbiting body. The 6 classical orbital elements were 
chosen for their physical significance in terms of how the satellite’s orbit is being 
affected by various perturbing forces. The particular satellite propagation theory 
being used here is loosely based on Brouwer's theory [7-11] for a future 
comparison basis with existing orbit improvement procedures. The detection 
scheme is an earth-fixed planar radar, the crossing of which causes power to be 
reflected back to one or more of multiple earth-based receivers [12], 
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The most poorly modeled phenomenon affecting low earth orbits is the 
problem of atmospheric drag. The overwhelming consensus [13-15] is that even 
the most complex atmospheric density models can only be consistently accurate to 
within about 20%. The Kalman filter incorporates a very simple decreasing 
exponential density model [16] corrupted by noise to track drag perturbations. 

The extended Kalman filter is used in many situations for nonlinear parameter 
estimation [17-19] and has been applied to the field of orbit estimation in the form 
of orbital parameter improvement by differential correction techniques [20,21], It 
is known that such a filter works well when observational data are plentiful 
[22,23]. However, proper application of such a filter to a highly nonlinear, 
minimal observation scenario, such as that posed by the fixed radar plane detection 
scheme, can also successfully "track" satellites. 
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III. TWO-BODY ORBIT FORMULATION 



Two-body orbital motion for a spherical earth can be described by a set of six 
parameters (three-dimensional, second order problem). Traditionally, these 
parameters take one of two forms. One method describes satellite motion in terms 
of Cartesian position and velocity. These vectors are propagated according to the 
two-body non-linear equations of motion. The other takes advantage of the fact 
that the equations of motion describe an orbital path which is a conic section, 
which, for earth-orbiting satellites, is always an ellipse. Satellite motion can then 
be analyzed as a mean angular displacement along that elliptical path. The second 
method is used here because of the linear nature of the mean angular motion and 
its clarity of physical significance. 

A. KEPLER TWO-BODY DYNAMICS 

The orbital path is described as an ellipse of fixed shape and size located in a 
plane of fixed orientation with respect to the earth's inertial reference frame. The 
quantities describing this path, the classical orbital elements [1], consist of five 

constants and one time-dependent quantity. Two of the constants (a, e) describe 
size and shape of the orbit, while the other three (i, Q, co) orient the orbital plane in 

space. The time-dependent quantity (v) pinpoints the actual location of the 

satellite on the orbital path. These elements are (see fig. 3.1): 
a (semi-major axis) - defines size of the elliptical orbit, 
e (eccentricity) - defines shape (non circularity) of orbit. 
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Figure 3.1. Classical orbital elements. 
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i (inclination) - defines the angle between equatorial plane normal (K) and 
orbital plane normal (W) . 

Q (longitude of the ascending node) - defines the angle between the vernal 
equinox direction (I) and the point where the satellite crosses 
the equatorial plane from southern to northern hemisphere 
(ascending node). 

oo (argument of perigee) - defines the angle between the ascending node and 
the satellite's closest point of approach (perigee). 

v (true anomaly) - defines the angle between perigee and the satellite's actual 
position. 

The elliptical orbital path, is described by the equation of a conic section 



1 + e cos v 

where r is the magnitude of the satellite's position vector at a specified v. In 
general, v is nonlinear with respect to time, so the mean anomaly M, which varies 

linearly with respect to time for purely Keplerian motion, will be used as a filter 
state instead. M is related to v through a quantity called eccentric anomaly E. It is 

defined as follows (see fig. 3.2): 



M = E - e sin E 



( 3 . 2 ) 



cosE - e 

cos v = 

1 - ecosE 



( 3 . 3 ) 



sin v = 



sin EVl -e 2 
1 - e cos E 



( 3 . 4 ) 
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Eqs. 3. 1-3.4 define the relationship between M and r . 

The filter state vector will be 

X k = [ 3 k e k *k ^k ^k] 

where 



' a k., ' 




a, 


e l-l 






ik + i 












W k+1 










a/h'.+Mj 



(3.5) 



(3.6) 



and Tk is the time between observations. For the classical 2-body case, it can be 
seen from eq. 3.6 that state propagation is a linear problem. 

In order to improve the estimate of the state vector using available 
observations, it becomes necessary to compare the states with the observables in a 
common reference frame. Satellite position, which is implicitly a function of the 
states through the relationship in eqs. 3. 1-3.4, can be written 



r 



p 



r cos v 
r sin v 
0 



(3.7) 



where r p denotes a position vector in the orbital (PQW) reference frame (see fig. 
3.1). The remaining states are used to transform satellite position from orbital to 
inertial (IJK) coordinates using 
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Figure 3.2. Eccentric vs. true anomaly. 
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where c0 = cos0 and s0 = sin 0 for brevity. (See App. A for development of 
C y *.) 
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Then satellite position in inertial coordinates will be 



r 1 =C 



Vr 



r cos v 
r sin v 



0 



(3.9) 



Now, inertial satellite coordinates need to be related to a measurement in the site- 
based reference frame. 



B. OBSERVATIONS 

The detection scheme modeled here is a fixed planar radar beam generated by 
multiple continuous wave transmitters along a great circle path (see fig. 3.3). As a 
satellite penetrates this beam, phase information is collected at one or more of the 
receivers. Antenna geometry at each site makes it possible to compute direction 
cosines, which are the pseudo-observations available for element set improvement. 
These observables are (see fig. 3.4) 



cos a = 



cos (3 = 



pJ3 

IpI 



p-N 

IpI 



(3.10) 



where East-West cosine = cos a and North-South cosine = cos p . (y can be 
calculated from: cos 2 y + cos 2 [3 + cos 2 a = 1.) 
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Figure 3.3. Radar plane geometry. 



We have 



pcosa = p-E 



p cos p = p • N 



(3.11) 



and 



cos a 
cos (5 



(3.12) 
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Figure 3.4. Receiver site geometry. 



is the vector of observables to be used for orbital parameter improvement. (See 
App. B for development of z in terms of the orbital elements.) 

The position vector in the site-based frame is 



r H =CV 



(3.13) 



where 



C H / = 



c(c§ 

-st{>c(az) - sf:c<f>s(az) 
s(az)s4>-c(az)s£c4> 



c£s<j) s( 

c(az)c4>-s(az)s£st() s(az)cf 
-s(az)c<t>-c(az)s£s(|) c(az)c£ 



( 3 . 14 ) 



and 

t — latitude of receiver 

4> = co e (t-t 0 )-lon nlf 
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CO*. = earth's rotation rate 

t 0 = time of last Greenwich crossing of vernal equinox 
lon site = longitude of receiver 
(See fig. 3.5) 

; App. A for development of C'-' 1 ). 

Recapitulating. 



both of which are, in general, nonlinear functions of the orbital elements \ t . 
Assuming for now that perturbations to the two body dynamics and observation 
noise can be modeled as additive zero mean, white noise processes, the state and 
measurement equations become 





(3.15) 



x k =f(x k . 1 ,T t _ I ) + w k _ 1 



(3.16) 



z k =g(* k ,T )+v k 



Kalman filter equations are then 





G k =P k/ H T [HP V H t + R 

* /k-l L /k-l 
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Xk X 



= x 



k X-i 





P ‘X=('- G .H) P «-, 



( 3 . 17 ) 
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where 



Q = E[w k w‘] 



R = E[v k vI] 



P 



k 




- x t 




T " 



0 = 






dx. 



‘K r T >-' 



H = 



5g(x k ,T k ) 



dx. 










(3.18) 






pcosa 

pcos(3J k 



(See App. B for derivations of <t> and H.) 

In words, the actual observations, z k , are compared with computed 
observations, z ^ j (based on the classical two-body propagation model), resulting 

in the Kalman filter innovations. These are multiplied by the filter gain, which is 
based on the dynamic model and measurement uncertainties and used to produce 
an improved estimate of the orbital elements. 
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C. RADAR PLANE INTERSECTION 



Meaningful application of the observations to orbital element set improvement 
requires the ability to predict subsequent penetrations of the radar plane by the 
satellite in question. A simple way to do this is to define the satellite's position 



out of plane component (Z) of position can be checked for a value of zero, 
yielding an in-plane condition. 

Satellite position, readily available in orbital frame coordinates (r p ) is 
transformed to radar plane coordinates by performing two coordinate 
transformations: one from orbital to inertial frame, the other from inertial to radar 
plane frame, or 



where r x denotes a position vector in the radar plane (XYZj reference frame and 



vector in a coordinate frame referenced to the radar plane (See fig. 3.6). Then the 




(3.19) 



C0C0 C0S0 s©^ 

rp rp rp 



= -s0 c0 0 

— S0 C0 -S0 S0 C0 

L rp rp rp 



(3.20) 



where 



0 = (O*(t-t o )-lon ll , 

©jpis the inclination of the radar plane to the equatorial plane 
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Ion, is the longitude where the maximum latitude of the radar plane occurs, 
which describes the location of the X-axis of the radar plane 
reference frame. (See fig. 3.6). 

(See App. A for development of C^.) 




Since the only component of satellite position of interest in determining radar 
plane intersection is the Z-component, the expression can be further simplified to 
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(3.21) 



r z = [0 0 l]C x ^C^r p 

= [-s0 ip P 1 ce-se ip P 4 s0 + ce ip P 7 ]rcosv 

+ [" se ^ p 2 ce “ se^PsSB + cG^P, ]r sin v 



where 



= 



P P 

r l r 2 

P P 

r 4 r 5 

P P 

r 7 r 8 



P 3 

P 6 

P 9 



(3.22) 



The in-plane condition of relevance to the filter is that which falls inside the field 
of view of one of multiple receiver sites, also placed along the great circle path 
circumscribed by the transmitters. The time epoch of this condition is the time 
close to which we expect to observe the satellite from one or more of the radar 
plane receivers. 
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IV. DOMINANT PERTURBATIONS 



Satellite orbital parameters deviate from Keplerian motion deterministically 
and randomly due to many factors. Ignoring for now such intentional human- 
caused phenomenon as plane change and stationkeeping maneuvers, the largest 
natural contributor to parameter variations is the nonuniformity of the gravitational 
field due to earth's oblateness. For low* earth orbits, the next most significant 
factor is atmospheric drag. 

A. GRAVITY PERTURBATIONS 

The largest contributor to perturbations to the Keplerian orbital elements 
arises from the oblateness of the earth. These perturbations take the form of a 
combination of secular and periodic variations of the classical orbital elements. 
Many solutions have been developed for this problem, each with its own strengths 
and weaknesses. All the solutions, however, are based on the gravitational 
potential 




where 



H* is earth's gravitational parameter 
r is the distance of the satellite from earth's center 
J n are n* order zonal harmonics 
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P n (sin 5) are Legendre polynomials 
5 is satellite declination 

and are more or less complex depending on the order to which the terms of the 
expression are carried. This expression for U already carries with it a degree of 
simplification in that it assumes axial symmetry for the earth. In truth, the equator 
is slightly elliptical, but the main effect is on the stability of geosynchronous 

satellites, since their orbits are in the equatorial plane [3]. Further simplification is 
accomplished by dropping terms for n >3. This is justified for medium accuracy 

applications because of relative magnitudes of the J n terms 



J 2 = 1.0826- 10 -3 



J 3 - —2 TO -6 



(4.2) 



J 5 = -2 TO" 7 



where J n successively decrease. So the simplified potential function is 




(4.3) 



Defining R = J 2 (3sin 2 5-l) as the perturbing function, the variation of that 
function with respect to each of the orbital elements is [ 4 ] 
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3R 

3a 

3R 

3e 




3J. 



2(1- -e 2 > 



■ [cos co(l + e cos v)(l -3sin 2 i ■ sin 2 u) 



-sin 2 i • sin 2u • sin v(2 + e cos v)j 

3R -3J, . . . . , 

— = — — sin i • cos i • sin " u 
3i r 



3R 

3co 




i - sin 2u 



3R 



= 0 



3D 
3R _ 1 dR 
3M ” n"dT 



(4.4) 



To develop a first order theory, the right-hand members of eq. 4.4 are 
expanded in a Taylor series about the initial orbital element values, retaining only 

the first term of the expansion, yielding expressions for the derivative of each 
orbital element with respect to true anomaly v. Integrating these expressions 

yields an analytical first order solution in the form 

x = x 0 +8,x(v) + 5 p x(v) (4.5) 



where 

x = instantaneous orbital elements 
x 0 = initial orbital elements 

8 s x = secular variations as a function of v 
5pX = periodic variations as a function of v 
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All six orbital elements experience periodic variations. However, only Q, to 



and M experience secular perturbations. These are 
5 s Q = -^[cosi 0 ]2(v-v o ) 



2 3J, 

S $ C0 = T-V 
2p‘ 



2--sin‘ i 
i 0 



( V “ V o) 



(4.6) 



5 M = n 



J, (< 



V 



1 % 



l+-r — 

a ; v r o j 



(l -3sin 2 6 o ) 



(t-t.) 



Note that is expressed in terms of independent variable t. which is the time at 
which true anomaly is equal to the value of v in 8 S Q and 5 s co. The periodic 
variations in to and M do carry 1/e 0 terms which approach singularity for near 
circular orbits. This can be explained by the fact that to and M are undefined for 

circular orbits and ill-defined for near circular orbits. Fortunately, the combination 
of co + M obviates the singular condition by effectively forming one new orbital 

element out of co and M. 

Long term propagation of the orbital elements, assuming a drag-free 
environment, is accomplished by applying only the secular variations, since 
periodic variations will have no net effect on the elements over time. However, 
the periodic variations will be seen later to come into play when observations are 
assimilated into the orbital element improvement process. Inclusion of these 
secular variations in the dynamic model yields 
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x k ^ = f(x k ,T k ) 



(4.7) 



u here 



a k 



e k 

K 



f( x k • T k ) = 



3J, 



co. + 



Q > -ttW'JIK*. - v k ) 

2pk 

3J, 



x ~ 2 

2p k 



-5. 2 . 

2 — sin i, 



K*i -v k ) 



(4.8) 



M k + n k 



J, 

l+-f 


(V) 




r . j 



>< 



(l -3 sin 2 5 k ) 



T, 



where 

Pk = a k (1 -e k ) is the semi-latus rectum of the orbit 

J 2 is the second order zonal harmonic (equatorial bulge). 

The gravitational force model which gives rise to these secular variations is 
relatively simple. It assumes that the earth is symmetric about its spin axis and 
about the equatorial plane. Noteworthily, the first three mean orbital elements 
experience no secular variation due to earth’s oblateness. 
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It may be possible to maintain identification of some satellites using secular 
variations alone. Continuous identification of a satellite, however, may depend 
upon the application of the periodic variations to the mean elements at observation 
times. These periodic variations are assimilated as an additive (superposition 
assumed) adjunct to the existing dynamic model and therefore are not technically a 
part of the propagation dynamics. The filter innovations resulting from a 
comparison of observables with predicted states are then applied, via the Kalman 
gains, to the mean elements, implying a prior removal of the periodic variations. 
These periodic variations.to first order, and ignoring terms of the order of e 0 and 
smaller, are 





+ — sin' i o -cos(2(0 o +3v 



12 






( 4 . 9 ) 
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B. DRAG PERTURBATIONS 



Satellites traveling in the altitude regime below 600 km experience non- 
negligible drag force due to atmospheric resistance [3]. This force is well 
represented as 



where 

m=mass of satellite 

C d =2 is a dimensionless drag coefficient 
A= effective cross sectional area of satellite 
V= velocity of satellite with respect to atmosphere 

p= local atmospheric density. 

1. Atmospheric Drag Effects 

To first order, drag force doesn't affect i and Q [6], The main effects of 

this drag force are to bring about secular decreases in a and e, i.e.. to decrease the 
size and eccentricity of the orbit. (O and M also experience changes due to drag, 

but these are negligible compared to the oblateness effects. 

It turns out that 




(4.10) 



m 



= —pa 

dM m Vl-ecosEy 



da _ C d A 2 /T + e cos 






(4.11) 



where 



24 



d(l - e cos E) 

T = 

1 + ecosE 

d = (l - e : Y 2 cos i 
n 



n = 




(4.12) 



and 



de 

dM 



Cn A /. 2 x (1 + ecosE)^ 2 

pa(l-e ) ^-(1-t) cosE 



m 



C D A . 
— - — paed 
2m 



(1 -ecosE)^ 
(\ -ecosE^ 



(4.13) 



1 1 + e cos E , 



(1 - T)sin 2 E 



All quantities in these two expressions are well known (according to the simplified 
model) except for atmospheric density p, the value of which is related to space and 

time in a complex manner. Atmospheric models delivering as good as 15% 
standard deviation are in the form of tables or very complicated empirical 

formulae or both [13]. For simplicity, it is advantageous to use an analytical 
function to model p. It is assumed that variations in p in the upper atmosphere 

(150 km to 750 km) are a superposition of four individual variations [4], 

The first variation is diurnal, or, daily. This fluctuation has to do with 
change in solar flux intensity associated with moving from darkness to daylight, 
and vice versa. Density, depending on latitude and altitude, may change by up to a 
factor of 10. 
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The second type is related to the 1 1 year solar cycle which is motivated 
by intensity of sunspot activity. Maximum density, occurring at the peak of 
sunspot activity, is roughly ten times the mean density. 

A third variation is semiannual, with density falling to a minimum in July 
and a maximum in October, with a less pronounced minimum and maximum in 
January and April, respectively. The range of this variation is, generally, less than 
1/20 the range of the 1 1 year solar cycle. 

The fourth type of variation is nonperiodic and unpredictable and is tied 
to geomagnetic storms rising from solar flare eruptions. Atmospheric density may 
increase by a factor of 10 from the geomagnetic quiescent value. 

Over the short term (days), barring magnetic storms, the diurnal variation 
dominates. However, over several years, the solar cycle has a more profound 
effect on atmospheric density. 

The modeling approach followed here w'as to start with a very simple, 
even crude, atmospheric density model, including the above considerations as 
necessary' to achieve desired accuracy. 

2. Simplified Drag Effects 

First, simpler expressions for changes in a and e due to drag can be 
developed by viewing these changes as a result of orbital specific energy being 
extracted by drag force [14]. Specific energy of an orbit, which is conserved for 
an unperturbed Keplerian orbit, can be expressed as 
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The change in a for an incremental change in c is then 



2a 2 

da = de 



(4.15) 



If it is assumed that most of the drag occurs near perigee (atmospheric density 
exponentially decreasing with altitude), then 



r = r = a(l -e) 



(4.16) 



from eq. 3.1, and 



V ' r = T-a, 
= UO + e) 



(4.17) 



from eqs. 4.13 and 4.15. For a satellite moving through dv near perigee (see fig. 
4.1), the specific energy extracted will be 




Figure 4.1. Satellite moving through perigee. 
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(4.18) 



dc = 



F 

— — rdv 



m 

2 



(C d A\ 



V m j 



pV : rdv 



Substituting eq. 4.17 into eq. 4.18 yields 



1 



de = — 

2 v m 



fC n A\ , 

'pn(l + e)dv 



(4.19) 



Substituting eq. 4.19 into eq. 4.15 gives an expression for decrease in semi-major 
a^is due to a decrease in orbital energy 



da = 




(4.20) 



Since the new orbit will pass through the same perigee point, all o f da shows up as 
a decrease in apogee radius, and 

dr p = da(l -e) - ade = 0 (4.21) 

by differentiating eq. 4.16. Solving for de gives 

de = — (1 -e) = - — p . A \(i -e : )pdv (4.22) 

a ^ m ) 



3. Atmospheric Drag Characterization 

A crude model for variation of p starts with the assumption that 

temperature and chemical composition of a gas remain constant [14], Then, 
density is directly proportional to pressure according to 
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p = kP 



(4.23) 



(Realistically, k is a function of temperature and chemical composition.) At an 
altitude h, the decrease in pressure accompanying an increase in altitude is 

dP = -pgdh (4.24) 

where g is acceleration of gravity. Since k is assumed constant. 

dP = -jjdp (4.25) 



Combining eqs 4.24 and 4.25 yields 

dp 

P 



= -kgdh 




(4.26) 



P = 



Po e 



-kgh 



So, in general, atmospheric density decreases exponentially with altitude. 
Eq. 4.26 can be rewritten 

P = P„e% (4.27) 

where p 0 (reference air density) and h 0 (density scale height) change for different 
altitude regimes. For example, using standard fixed values of p from the "U.S. 



29 



Standard Atmosphere, 1976," a reasonable fit can be made by applying the values 
in Table 4.1 for the ranges of altitude shown. 



Table 4.1. DRAG PARAMETERS. 



Altitude 


p 0 (kg/m 3 ) 


zo (km) 


80 km<z<120 km 


18.5 


5.8 


120 km<z<140 km 


1 .2 1 x 10' 3 


11 


1 40 km<z<200 km 


3.00 x 10-6 


21 
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V. THE RECURSIVE ORBITAL PARAMETER ESTIMATOR 



An extended Kalman filter, designed for orbital parameter estimation, is based 
on the earth-fixed radar plane detection scenario described in Chap. Ill (Also see 
fig. 5.1). The filter propagates the orbital elements in time according to a first 
order nonlinear dynamic model which takes into account the dominant 
perturbations arising from earth oblateness and atmospheric drag. Upon receiving 
observations, in the form of pairs of direction cosines, the filter calculates an 
improvement to the current parameter estimate. The filter operates recursively as 
observations become available. The filter is written in Matlab code and can be 
found in App. C. 

A. FILTER INITIALIZATION 

The estimator is designed to receive two data files. One contains an initial set 
of orbital elements along with subsequent pairs of available direction cosines (see 
App. C). The other holds information about the receiver sites. Values of pertinent 
constants are also set, as are the initial times necessary for orbit propagation and 
estimation of observables. 
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Figure 5.1. Estimator flow chart. 
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1. Input Data Files 

Two input files are required for filter operation. One contains an initial 
set of orbital parameters, along with subsequent sets of direction cosines. This 
information is all tagged with corresponding times (see App. C). The last five 
orbital elements occur in classical orbital element form, but, in lieu of semi-major 
axis, orbital period is given. The conversion is simply 



a = 



271 
V 

(T 

sat 

v 271 J 



(5.1) 



w here T sa[ is the period of the satellite. Note here the apparent disappearance of 
earth's gravitational parameter This convenience is brought about by the use 



of canonical units [1], where 



i 3 



=1 



DU\ 

tuT 



1 DL T £= 6378.135 km and is earth's mean equatorial radius. 

1 TU e = 13.44686 min and is the time it takes for a satellite traveling in a 
circular orbit at a radius of 1 DU to travel 1 DU of arc length. 

Canonical units are used in all filter calculations. The other input file contains 

information on the receiver site locations,altitudes, and azimuths. The azimuthal 

measurement indicates, for each receiver, the angular offset of its coordinate frame 

from true north. 
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2. Constants 



The earth's apparent rotation rate is initialized as 

(o* = .05883378 (5.2) 

The radar plane’s inclination to the equatorial plane and the longitude of its center 
line are, respectively, 

6^ = 33.58310° 

(5.3) 

Ion, = -101.31348° 

The radar plane is offset from earth's center by about 20 km, or 

r.,,.., =.0031 DU. (5.4) 



3. Time 

All times are given in Universal Coordinated Time (UCT) and converted 
to TUeS for filter use. Tg rnv , C h,which is time elapsed since the Greenwich meridian 
last passed the geocentric inertial I axis, is computed from a reference angle 
(measured from I to Greenwich) at 0 hrs, 1 Jan 92, of 

0 c =98.920250931“ (5.5) 

gnwcn v x 

using the co^. given previously. 
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4. Filter Matrices 



The error covariance matrix P is initialized to reflect the level of 
confidence placed in the accuracy of the initial orbital element set. The 
measurement error covariance matrix R is set to 

R = [°T J } (5.6) 

L 0 

where o c0fa = 0 ^ = .0002° reflects the uncertainty of the available direction 
cosines. 

B. STATE PROPAGATION 

The filter states are propagated according to the two-bodv Keplerian orbital 
dynamics perturbed by earth oblateness and atmospheric drag. The nonlinear state 
equations are 







2 — sin 
2 



-Vi 

(l -3 sin 2 8 k ) T k 




(5.7) 
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and are propagated at arbitrary increments of • This is continued until a 



crossing of the radar plane is detected by a sign change of the radar-plane-normal 
component of satellite position. 

C. RADAR PLANE INTERSECTION 

Detection of a radar plane crossing calls a function which iterates to an in- 
plane condition for estimated satellite position. This is accomplished using a 
Newton-Rapson iteration on the expression for the plane-normal component of 
satellite position. The iteration calls for a calculation of an approximate AM 
which will drive r z to near zero. This is given by 



A 




(5-8) 



where, from eq. 3.21, 



t z = X , r cos v + X , r sin v 



(5.9) 



and 





3M 0M 



(5.10) 




OM 



Taking the individual derivatives yields 
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ax, 

avT 

ax 2 

~3m 



= -P. sin 0 T 



= -P„ sin 0^ 

2 rp 



a(cos 0 ) 

0M 

a(cos 0 ) 

aM 



- P, sin 0 m 

4 ip 

- P 5 sin 0 

5 rp 



3(sin 0) 

a.\i 

a(sin 0) 

aM 



where 



and 



where 



a(cos 0) 

0M 

3(sin 0) 

a\i 



= -co 4 .a"^(sin 0) 
= co e a"^(cos 0) 



3(r cos v) 


a sin E 


a.\i 


1 -e cosE 


3(r sin v) 


aVl - e : cos E 


3M 


1-ecosE 



rcos v = a(cosE-e) 
r sin v = aVl -e : sin E 

aE _ i 

3M 1 - e cos E 

This iteration is found to converge to within z r <£ in three to 
w'here e = 10"\ 



( 5 . 11 ) 



( 5 . 12 ) 



( 5 . 13 ) 



( 5 . 14 ) 



four iterations, 
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The function then determines whether the satellite is within the field of \iew 
of one or more receiver sites by calculating the in-plane r x and r y components of 
satellite position and checking the inequality (see fig. 5.2) 

—=====: > cos 20' (5.15) 

V r » +r ; 

It should be noted here that the field of view needed to be opened up to pick up 
some satellite observations by receivers at the edges of the radar plane. 




Figure 5.2. Sensor field of view. 



D. ASSIMILATION OF OBSERVATIONS 

Once a satellite is in the field of view of radar plane coverage, the input data 
file is checked for observations at or near its time of arrival. If a pair of direction 
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cosines is available, then a measurement estimate is calculated according to the 
development in App. B, where 



z = 




( 5 . 16 ) 



<t> is also calculated according to App. B. The plant noise covariance matrix Q is 
calculated using the squares of expected levels of uncertainty in each of the states. 
For example, if the two-body problem is the dynamic model, then the Q matrix 
would be 

~c; 0 0 0 0 

0 o 2 0 0 0 

0 0 o 2 0 0 

Q = , 

0 0 0 o 2 n 0 

0 0 0 0 o 2 

0 0 0 0 0 

where c n are proportional to the magnitude of the largest unmodeled perturbation, 
in the two-body case, earth oblateness. So, for the two-body problem, 



0 

0 

0 

0 

0 



( 5 . 17 ) 
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3J ; a 



(l - e 7 ) sin : i • cos[2(co + v)] 



3J, (. . 



1 . 



o. 

Or 



-V 



1 - — sin ‘ i cos v + — sin 2 i • cos[2co + v] 



3J : 

8 ? 



sin 2i • cos[2(co + v)] 



3J, 



-V 

3J, ( 



r- n cos i ■ T. 



2 P 
3J, ( 



2p‘ 



2 --sin' i T k 



- ^sin : ij(l -e : )^T. 



(5.18) 



Now. the error covariance matrix can be calculated as 

I\ /k _. =4)P«. X O t +Q (5.19) 

H is then calculated as given in App. B, making possible the calculation of filter 
gain 



G k =P^_ j H t (hP^_ i H t +r)’ S (5.20) 

which is then applied to the residuals z-z to calculate improved orbital 
parameters \y 

It should be noted that, throughout the filter, the quantities v and E are needed 
for various calculations. E is calculated fom M using Newton-Rapson iteration on 

M = E - e sin E (5.21) 



and v is calculated by solving for 
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( 5 . 22 ) 



v = sin 



\1 - e 2 sin E 
1 - e cos E 



then using 



v = cos 



cosE-e " 
1-ecosE. 



( 5 . 23 ) 



to resolve the ambiguity. 
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VI. RESULTS 



This extended Kalman filter is designed to estimate the orbital parameters of a 
satellite, given observables in the form of a pair of direction cosines obtained upon 
the satellite's passage of an earth-fixed, stationary planar radar beam. One day's 
worth of data for each of three satellites was obtained for filter testing and 
adjustment. Results of the application of three successively more comprehensive 
filter models are compared with each other. 

The first filter models only two-body motion, the second adds secular 
perturbations due to earth oblateness, and the third also includes periodic 
perturbations due to oblateness as well as atmospheric drag effects. Each filter, of 
course, models unknown perturbations as zero mean, white, Gaussian noise, in the 
tradition of Kalman filtering. Unfortunately, none of the available satellite data 
include an orbit in the drag regime. The filters will be compared in the areas of 
radar plane time of arrival, a residual small sample statistic, and orbital parameter 
estimation. 

A. RADAR PLANE TIME OF ARRIVAL 

With orbiting objects arriving at the radar plane on the order of every few 
seconds, it is important that prediction of arrival times be fairly accurate. Current 
methods yield accuracies of within a second. The best that can be hoped for under 
the current detection configuration is about 0.25 second, which is the sensor time 
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uncertainty. After some crude tuning of the three filters, a comparison of arrival 
times for the three satellites is shown in Table 6.1, where At = t obf -t B1 . 

It can be seen that, in all cases and for all models, the first detection time after 
the estimator is started up (which occurs on the order of 12 hours after filter 
initiation) is always mispredicted by the filter by a time on the order of seconds. 
However, it is seen that the secular and periodic models cut down that time error 
by a factor of from three to four. Subsequent detections yield differences in time 
between calculated and observed plane intersections of less than 1 second for a 
one-orbit-later detection and between two and three seconds if detection occurs 
next on the "backside” of the orbit (about a half day later). Although this is true 



Table 6.1. OBSERVED VS. ESTIMATED TIMES OF ARRIVAL. 
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for the two-bodv and secular cases, the incorporation of periodic perturbations 
yields subsequent time differences of less than one second in all cases. 

B. RESIDUAL SMALL SAMPLE STATISTIC 

The filter residual is the difference between actual and estimated observations, 
and, as such, is the basis for a good performance measure. The small sample 
variance of the residual is calculated for each case as follows: 

o^X — -X ! (6.1) 

n 

where 

X = I- (6.2) 

i-i n 

is the small sample mean, x, being the i* residual. Plots of these variances for each 
case are shown in figs. 6.1 - 6.3. 

It is evident that the residuals for the N-S direction cosines are always very 
small. This is because this direction cosine is always close to 90° and the 
geometry of the problem always forces the estimate of this direction cosine to be 
very close to 90°. However, the E-W (in-plane) cosines are more interesting to 
observe. In general, the trend is toward smaller residual variances from 
observation to observation and from simple to more complex dynamic model. The 
only case where this does not hold true is for Satellite #117, where it can be seen 
that the two-body model is sufficiently lost after approximately eight hours to 
yield higher residuals than at the previous set of observations. Note that a "stack" 
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of asterisks denotes several observations gathered at the same radar plane crossing, 
which are then applied one at a time for parameter improvement. 

It can be concluded from this limited residual data that, for the relatively short 
times included in the data, the filter is performing acceptably. 

C. ORBITAL PARAMETER ESTIMATION 

The ultimate objective of this filter is to be able to predict what a satellite's 
position will be at some future time. Hopefully, the filter is predicting plane 
intersection times by accurately estimating orbital parameters. Figs. 6.4 - 6.6 
show time propagating values of the four parameters that are expected to be 
varying most with respect to their two-body propagation values. All three 
satellites show similar results. In all cases, semi-major axis estimates differ 
greatly from model to model. Eccentricity estimates for the two-body and secular 
models are very similar, with the periodic model showing the greatest deviation. 
All estimates are very similar for the longitude of the ascending node, though, 
between observations, the two-body model does not have a mechanism to 
accurately propagate this element. Argument of perigee estimates are very close 
between the secular and periodic models, while the two-body makes very little 
correction to this parameter. 
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Figure 6.1. Variance of residuals (sat. # 116). 
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Figure 6,2. Variance of residuals (sat. # 117). 
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Figure 6.3. Variance of residuals (sat. # 118). 
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Figure 6.4. Orbital element estimates (sat. # 116). 
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Figure 6.5. Orbital element estimates (sat. # 117). 
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Figure 6.6. Orbital element estimates (sat. # 118). 
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VII. CONCLUSIONS AND RECOMMENDATIONS 



Orbital parameter estimation was attempted using an extended Kalman filter 
for each of three orbital dynamic models. Observations were angles only and 
made roughly twice a day per satellite passing through the radar plane. It was 
found that the two-body orbital model, including secular and periodic variations 
due to earth oblateness, was able to predict subsequent radar plane crossings to 
within a second of actual crossings. 

The filter, which includes a crude model for atmospheric drag effects, was not 
tested on any satellites in the drag regime. Future research should include data 
from satellites with altitudes less than 600 km. Also, longer term data should be 
obtained to check filter performance over the long haul and also to see if the filter 
could track longer term disturbances or perturbation; to the orbit under 
consideration. 

Tuning of the filter could also be accomplished by obtaining a larger data base 
of satellite orbital elements and corresponding observations. This could be done 
manually, but it would be more exciting to see if an artificial neural network could 
be designed and trained to accomplish the same task. 
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APPENDIX A 



COORDINATE TRANSFORMATIONS 

The calculations and algorithms presented in this paper necessitate describing 
vectors in several different coordinate frames. These frames are all orthogonal. 
Transformation of a vector from one reference frame’s coordinates to another's is 
accomplished using direction cosine matrices (DCMs). 



Reference Frames 

The reference frames of interest are as follows (see fig. A.l): 

Geocentric inertial (IJK) - The I vector points inertially in the vernal equinox 
direction, with I and J lying in the equatorial plane, and K 
piercing the north pole (coincident with earth's spin axis). 

Orbit-fixed (PQW) - The orbital plane is the fundamental plane, with P 
pointing to perigee (point of closest approach), Q is the semi- 
latus rectum direction, and W is the orbit-normal (also orbital 
angular momentum direction). 

Radar plane (XYZ) - The radar plane is the fundamental plane, with X 
pointing to the maximum north latitude point, Y lies in the 
equatorial plane, and Z is the radar plane normal. 

Site-based (HEN) - H is local vertical, E is east, and N is north. The origin 
coincides with the radar site location of interest. 



Direction Cosine Matrices 

Three DCMs are developed here, performing one orthogonal rotation at a 
time, then combining the rotations. 
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UK (Geocentric inertial) 




PQW (Orbit-fixed) 





Figure A.l. Reference frames. 



Orbit -fixed to Inertial 

Satellite position is most simply expressed in orbit-fixed coordinates, 
but for comparison with earth-fixed quantities, it is useful to transform both to an 



56 



intermediate reference frame. The inertial coordinate frame (in this case, 
geocentric inertial) is the obvious choice. Fig. A. 2 shows the individual rotations 
required to transform a vector from orbit-fixed to inertial coordinates. 




Figure A.2. Orbit-fixed to inertial rotations. 
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Combining the rotations yields 
= [rot3][rot2][rotl] 



ccocQ - scosQci 


ccosQ + scocQci 


scosi 


-scocQ - ccosQci 


-scosQ + ccocQci 


ccosi 


sQsi 


-sicQ 


ci 



Because the transformation is orthogonal, 

c l/p =[c p ^] T = 



ccocQ - scosQci 


-scocQ - ccosQci 


sQsi 


ccosQ + scocQci 


-scosQ + ccocQci 


-sicQ 


scosi 


ccosi 


ci 



P P 

r, J 3 



P, P 5 

P, P 8 



(A.l) 



(A. 2) 



Inertial to Site-based 

Fig. A. 3 shows the individual rotations required to transform a 
vector's coordinates from the inertial to a site-based reference frame. Rotation #1 
keeps track of earth rotation. Rotation #3 points out the fact that the E and N 
vectors are not truly east- and north-pointing. E is in fact tangent to the great 
circle circumscribed by the radar plane, and N is orthogonal to E in the local 
horizontal plane. 
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where (?= co$(t — t c )+ lon cti 

is earth's rotation rate 

t c is time of last Greenwich crossing of I axis 
lon tlU . is longitude of radar site 



K 




rot 2 = 



cos( 


0 


sin i 


0 


1 


0 


-sin L 


0 


cosl 



where Y is latitude of radar site 




Figure A.3. Inertial to site-based rotations. 
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Combining the rotations yields 



C^ ; = [rot3][rot2][rotl] 



(A. 3) 



c(c§ 



c£s <t> 



-s<{)c(az)-s^c4)s(az) c4>c(az) - s(az)s/sb s(az)c£ 
s<{>s(az)-s£c4)c(az) -c<{)s(az)-s£s<t>c(az) c(az)c£ 



Inertial to Radar Plane 

The transformation from inertial to radar plane coordinates is very 
useful to the determination of satellite crossings of the radar plane. Rotation #1 
accounts for earth rotation. Fig. A. 4 shows the individual rotations necessary to 
perform this transformation. 

Combining these rotations yields 



C0C0 c0 s0 s©^ 

rp rp rp 



C^ = — s0 c0 0 

~S0_C0 -S0 ro S0 C0_ 

L rp rp rp 



(A. 4) 
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where 0 = co®(t — 1 0 )+ Ion, 

Oe is earth's rotation rate 

t 6 is time of last Greenwich crossing of I axis 

Ion, is West longitude of maximum North 
latitude of radar plane 



K 




rot 2 = 



COS0 


0 


sin 0 


0 


1 


0 


-sinG* 


0 


COS0 



where 0^ is inclination of radar plane to equatorial plane 



Figure A.4. Inertial to radar plane rotations. 
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APPENDIX B 



LINEARIZED FILTER MATRICES 6 AND H 

The extended Kalman filter can be viewed as a predictor-corrector type 
estimator. States are predicted using non-linear plant dynamics. Then a correction 
is computed by applying the Kalman gain to the filter residuals (observations 
minus estimates). Calculation of the Kalman gain requires some type of 
linearization of the plant and measurement dynamics (see eq. 3.18). Normally a 
Taylor series expansion is used, keeping only first order terms which are evaluated 
at the best current state estimate. The general form of the linearization of a 
nonlinear discrete function 




(B.l) 



. f n,(x k ,T k )_ 



is 




df : df t 




(B.2) 
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Development of Linear Plant Matrix 

A. 

As will be seen, the state matrix <J> can be as complex as one lets it 

become. Following is a presentation of three progressively more complex models; 

1) Keplerian two-body motion, 

2) inclusion of dominant first order earth oblateness gravity effects, and 

3) dominant atmospheric drag effects. 

Two-Body Problem 

State dynamics for a Keplerian 2-body orbit are given as 



’a I+I ' 




a. 


e „, 






* k+1 




L 








w k+1 






.M„. 




A+a;^. 






Linearizing gives 



1 0 0 0 0 0 

0 1 0 0 0 0 

0 0 1 0 0 0 

0 0 0 1 0 0 

0 0 0 0 1 0 

-~a' k P T k 0 0 0 0 1 

2 k 



(B.3) 



(B.4) 
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Secular Variations Due to Earth’s Oblateness 



To first order, the last three orbital elements are affected by secular 
variations, so that 



where 



x ... = 



3J 



Q k - — 4(cosi k )Av 



2 Pl 



CO, 



M k +a 7 S 



3J. 



2p‘ V 



(, 5 . 2 . 
2 — sin 



Av 



1 + - : , * (l — 3 sin 2 5 k ) 



% 



T, 






(B.5) 



r k = 



l + e k cos v k 



sin§ k = sin i k sin(\' k +co k ) 



(B.6) 
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Linearizing gives 
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(B.7) 



where the terms in the fourth row are found to be 

an,., 3J 



3a, 


P k 3 


an.., 


_ 3J 2 


de k 


pI 




3J 2 


ai, 


"2 P 2 


an,., 


3J 2 


dM k 


~2 P 2 



, dp 

cosi k Av 



3a 

i 

dp i 



cosi k Av^ i - : - 

de. 



sin i k Av 



cost 



: jK. 

k dM, 



(B.8) 
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where 



dp* 



5a k 



= 1-e: 



dPk 



5e k 



-2ae 



5v k l(l + e 2 )sinE + 2esinEcosE 
5M k -sin v k (1 -e cos E) 3 



The fifth row terms are found to be 



dco k+i 

5a k 






5e k 





5co k+ , — 1 5 J 2 . . 

= sin i . cos i k Av 

5i k 2p t 



5co k „ ^ -3J 2 a 5v k 
5M k 2p k 3 * 5M k 



where 



A 
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i 



(B.9) 



(B.10) 



(B.l 1) 
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Sixth row terms can be evaluated by 



where 
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(B.12) 



(B.13) 



67 



and 



A . =1—3 sin 2 5 ^ 
A s# =l-e; 



A 6 

6 <t 

A 7 

' <t> 

A 8 



>2 



1 + "T A 4„ 



A 



= l + e k cos v k 



J, 



= -f (l-3sin 2 8 k ) 



(B.14) 



Development of z and H 

The relationship between the measurements (pseudo-observables) and 
orbital parameters (filter states) is contained in the nonlinear function 



= g( x k ’T k ) 



(B.15) 



and is the same regardless of the state dynamics model being considered. The 
observables of interest here are the north-south and east-west direction cosines, or 



cos a 
cos[3 



(B . 1 6 ) 



and are related to the range vector p by 







Pe / 


cos a 




/P 


cos(3 




Pn/ 






L /pJ 



(B.17) 
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Note (from fig. B.l) that the HEN coordinate system in which p is received is not 
truly up, east and north. H is local vertical, but E is tangent to the radar plane 
great circle, and N is orthogonal to E in the local horizontal plane. This site- 
dependent information is included in the transformation matrix (see App. A) 




So p E , p N , and range magnitude p need to be expressed in terms of the 
orbital parameters in order to form z = g(x, T). It can be stated that 





’p H " 




pcosy* 


p = 


Pe 
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pcosa 




_Pn_ 




pcosP 
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In order to describe the direction cosines in term of orbital parameters, first 
consider fig. B.2. Earth’s flattening has been exaggerated in the figure to clearly 
show the effect on radar plane orientation. From the figure, it is clear that 

r = r ofr«; + +P (B.19) 



where 



R*. = R(lat) + alt 
R(lat) = 

i 1 



1 __ 


2 ~ ^298.26 


1 


298.26 



sin : (lat) 



(B.20) 



r offs« = 20 km 




Figure B.2. Relationship between satellite position and range. 
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R(lat) is the latitude dependent radius of a point on the reference geoid (due to 
flattening) [7] and is assumed to be in the H direction. r 0 ff se tis assumed to be 
totally in the negative N direction, and alt^e is the receiver’s altitude. Combining 
the elements of eqs. B.18 and B.19 yields 



r 



H 



R au + pcos y 
pcosa 

“ r offs« +pcosp 



(B.21) 



r H is the transformation of r p from the orbit-fixed coordinate system to site-based 
coordinates, or 



= C^C^ p r p 



H, H 2 H 3 

h 4 h 5 h 6 
_H 7 H s h 9 J[p 7 P g 

r(R 7 cos v + R s sin v) 
r(Rj cos v + R 2 sin v) 
r(R 4 cosv + R 5 sin v) 



r cos v 
r sin v 
0 



(B.22) 



where 

R,=H,P, +H S P <+ H,P, 

R 2 =H,P i+ H s P s +H 6 P, 

R, =H,P, +H,P 4 +H,P, 

R s =H,P 1 +H,P s +H,P, ( b - 23 ) 

R, =H,P, +H ! P,+H,P, 

R,=H,P ! +H ! P s + H 5 P, 
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Combining eqs. B.18,B.21, and B.22 yields 



"Ph" 




r(R 7 cos v + R s sin v) - R aic 


Pe 


= 


r(R t cos v + R, sin v) 


.Pn. 




_r(R 4 cos v+ R s sinv) + r offseI _ 



(B.24) 



From eq. B.24, p can be calculated as 



P = 



VP H P E P N 



(B.25) 



Now, eq. B.15 can be used as the estimate of the observation to compute filter 

residuals. But the Kalman gain depends on the matrix H, which is derived by 
linearizing g(x k , T k ) with respect to x K . 

The form of H is as follows 



a 





0a 




3M J 



(B.26) 
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where 
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Following are the details of this development. 
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(B.27) 



(B.28) 
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where 



(B.29) 



3r 



1-e' 



3a 1 + ecosv 



The expressions for the change in p with respect to eccentricity e are identical to 
eq. B.28 except for the replacement of with which is 



dr -a(2e + cos v + e 2 cos v) 
3e (1 + ecosv) 2 



(B.30) 



Since the R n terms in eq B.24 are functions of the orbital elements i, Q, and to, 
the expressions for the change in p with respect to these three elements are as 
follows, with the replacement of 3C by the appropriate element: 
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(B.31) 
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The expressions for the change in p with respect to a change in mean anomaly 



M is more involved because of the relationship between mean and true anomalies. 
This relationship is repeated here for convenience: 

M = E - e sin E (B.32) 



cosE - e 

cos v = 

1 - ecosE 



(B.33) 



sin v = 



sin E Vl - e 2 
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(B.34) 



Taking partial derivatives as appropriate yields 
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where 
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(B.36) 



All the expressions for the partial derivatives of the R n terms can be found in 
Table B.l. 
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Table B.l. DEVELOPMENT OF TERMS FOR H. 
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APPENDIX C 



FILTER CODE (MATLAB) 

Following is the Matlab code for the extended Kalman filter described in this 
paper. This particular filter models the 2-body orbital dynamics problem, 
including secular and periodic perturbations. It requires 2 data files for input (see 
header of filter code): 

datain - contains a time-tagged initial orbital element set and an unspecified 
number of time-tagged pairs of direction cosines and 
corresponding receiving sites 

siteinfo - contains pertinent information about the set of receiver sites 
comprising the radar plane detection system. 

The filter calls 3 functions: 

enewton - performs a Newton-Rapson iteration to obtain eccentric anomaly, 
given a mean anomaly and eccentricity 

intper - performs a Newton-Rapson iteration to propagate satellite position to 
an in-radar-plane condition. 

tH'opi - forces angular measurements to a number between 0 and 2k. 



Extended Kalman Estimator 



The code for these functions follows the filter code. 



This estimator requires 2 input 
DATAIN & SITEINFO. 



datain= [ 250 116 82 191626.933 

103.591322 .00756997 66.81264 

116.53101 

83 65024.876 6 14.839 .092 

83 83436.223 5 -76.912 .012 



0 0 

189.71773 

0 

0 



26 



files 

.77689 



to 



run 
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83 


83436. 


259 


4 


-63.554 .009 


0 


83 


83436. 


397 


1 


58.962 .062 


0 


83 


83436. 


406 


3 


-57.058 .043 


0 


83 


83436. 


425 


6 


-78.987 .034 


0 


83 


160057 . 


722 


6 


64.179 .030 


0 


83 


174304 . 


847 


3 


-16.315 .020 


0 


83 


174304 . 


904 


5 


-58.495 .007 


0 


83 


174304 . 


906 


6 


-62.516 -.001 


0 


83 


174304 . 


908 


1 


74.290 .012 


0] / 


sit 


einf o= [ 1 


32 . 


57840 -116.97016 


1 . 87473e-5 


2 


33.44600 


-10 


6 


.99824 2 . 21244 9e 


-4 3.13900 


3 


33.33239 


-93 


. 


55039 8 . 1 4 31 07e- 


6 -4.27563 


4 


33.14672 


-91 


. 


02096 6.738134e- 


7 -5.66895 


5 


32.28766 


-83 




53628 1. 11824 3e- 


5 -9.71924 


6 


32.04323 


-81 


• 


92609 3.866064e- 


6 -10.5800 


* * x 


*•*■*■*★*■★★ 


* * * * * 


★ 


**■*★*■*'***★'*'**■*■*'*'* 
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function [xkeep, Tkeep, Tint , del Tint , res id, xint , Pkeep] = . . . 

orbper (datain, siteinfo) / 
inpels= [data in ( 1 , : ) ' ; data in (2, : ) ' ] ; 
observ=dat a in (3 : length (data in ( : , 1 ) ) ,1:5) ; 

K=inpels ( 1 ) ; 
for i = l : 7 

lat (i)=siteinfo(i,2) * pi/1 80; 
lon(i) =siteinfo(i, 3) *pi/180; 
alt (i) =siteinfo(i, 4 ) ; 
az (i)=siteinfo (i, 5) *pi/180; 
end 

S'* , x , ***'*'*'rt******'****'*'*'*'*'*'*'******** , **'***'rt'***'*'***'*'*xx*'*7fr'*xx-*-*'rt'*'*'*** 

o 

THIS PROGRAM PROPAGATES A SATELLITE ORE IT AS A 2-BODY PROBLEM 
WITH SECULAR VARIATIONS IN OMEGA, omega, and MEAN ANOMALY. IT 
CALLS THE FOLLOWING FUNCTIONS: 

TWOPI - places angle between 0 & 2*pi. 

INTSEC - iterates to an in-fence-plane condition when a fence 
plane crossing is detected. 

ENEWTON - solves for eccentric anomaly, given values for mean 
anomaly and eccentricity. 



* DEFINE EARTH ROTATION RATE * 

earthrot= . 0588337 8171 654 / % Define earth rotation rate (rad/TU) . 

J2=1.082645e-3; 

* DEFINE VARIOUS TERMS * 

R=[4e-8 0/0 4e-8]; % Measurement noise matrix. 
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Pkkml=eye (6) *le-8; % Initialize error covariance rr.acrix. 

Pkeep ( : , 1 : 6) =Pkkml ; 
resid ( : , 1) =[0/ 0] ; 

G=zeros ( 6, 2 ) ; % Initialize Kalman gain. 

finc=33 . 58310*pi/l 80; % Define fence inclination. 

cosfinc=cos (fine) ; 
sinfinc=sin (fine) ; 

lonx=-101 . 31348*pi/180; % Define location of fence x-axis. 

* INITIALIZE STATES (0R3ITAL ELEMENTS) * 

xnow= [ ( inpels (7) / (2 *pi * 1 3 . 44 68 9317 ) ) A (2/3) ; 

inpels (8) ; inpels (9:12) .* pi/180 j ; 
xkeep ( : , 1 ) =xnow; % Save original states. 

* CALCULATE INITIAL ECCENTRIC ANOMALY * 

O 

EA ( 1 ) =Enewton (xnow ( 6) , xnow ( 6) , xnow (2 ) ) / 

o 

* FIND INITIAL r and TA * 

G 

cosTA= (cos ( E A ( 1 ) ) -xnow (2 ) ) / (l-xnow(2) ’cos (EA(1) ) ) ; 
sir.TA=sqrt ( 1 - xnow (2 ) A 2 ) ’sin (EA ( 1 ) ) / ( 1 - xnow (2 ) ’cos (EA ( 1 ) ) ) ; 
r=xnow ( 1 ) * ( 1 - xnow (2 ) A 2 ) / ( 1 i-xnow (2 ) * cosTA) ; 

TA ( 1 ) =twopi (asin (sinTA) ) ; 
if TA(l)<pi 
if cosTA<0 

TA ( 1 ) =pi -TA ( 1 ) ; % Perform quadrant check 

end % for true anomaly, 

else 

if cosTA<0 
TA (1 ) =3 *pi -TA ( 1 ) ; 
end 
end 

* CALCULATE TRANSFORMATION MATRIX (ORBITAL TO INERTIAL) 

cosperi=cos (xnow (5) ) ; 
sinperi = sin (xnow (5) ) ; 
cosascnd=cos (xnow (4 ) ) ; % 

sinascnd=sin (xnow ( 4 ) ) / % 

costilt=cos (xnow (3) ) / 
sintilt = sin (xnow (3) ) ; 

Pl=cosperi*cosascnd-sinperi*sinascnd*costilt; 

P2=-sinperi *cosascnd-cosperi *sina send* costilt / 

P4=cosperi * sinasend+sinperi’eosasend* cost ilt ; 



Calculate sines & cosines of 
angular orbital elements. 
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P5 = -sinperi * sinascnc+cosperi * cos as end* cost i It ; 
P7=sinperi * sincilt; 

P8=ccsperi*sintilt; 

5£-X*-K***'k*xx**'k*'k***’k'X****'**'**'*'*'*****'*'*****-**yr'*'*'rv r 'X'* : t * -r v 



* INITIALIZE TIME * 

S.*x****r********************* **************^***-jir-**-ir * r * * t 

c 

one Jan92=98 . 920250931 *pi/ 180; 

TUperday=l 07.0879334097483; 
secperTU=806 . 81359; 

JD=inpels (3) ; 
time=inpels (4 ) ; 
hrs = fix (time*le-4) ; 
min=fi>: (rem (time , le4) * le-2 ) ; 
sec=rem (time, le2) ; 

Tfilter= (hrs * 3 600+min* 60 + sec) / secperTU; 

Tgnwch=rem (earthrot *TUperday* (Tf ilter /TUperday+ . . . 

(JD-i) ) +oneJan92, 2*pi) /earthrot; 

Tcld=Tf ilter ; 

Tkeep (1) =JD+Tf ilter * secperTU/ 8 64 0 0 ; % Save original 



-rir 






. ime 



* INITIALIZE TIME OF FIRST OBSERVATION 



time=observ (1, 2) ; 

hrs=fix (time*le-4) ; 

min=f ix (rem (time, le4 ) * le-2 ) ; 

sec=rem (time, le2 ) ; 

timeobs= (hrs* 3 600+min* 60 + sec) /secperTU; 

* CALCULATE OUT-OF-FENCE-PLANE COORDINATE OF SATELLITE POSITION * 

O 

theta=twopi (earthrot *Tgnwchi-lonx) ; 
sintheta=sin (theta) ; 
costheta=cos (theta) ; 

zfence=.0031+ (-sinfinc*costheta*Pl . . . 

-sinf inc* sintheta*P4 . . . 

+cosfinc*P7) *r*cosTA. . . % Calculate current magnitude 

+ (-sinf inc*costheta*P2 . . . % of out-of-plane coordinate, 
-sinf inc*sintheta*P5 . . . 

+cosfinc*P8) *r*sinTA; 



* PROPAGATE ORBITAL ELEMENTS TO NEXT FENCE-PLANE INTERSECTION * 

J2=l . 082645e-3; 

P=l; 
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n=0 ; 

f cr k=l : K 

Tloop=pi/8*xnow (1) *1.5; 

^****X************rt**xx*XX*XXXX*X*xx*X*x*X**x**XxxxxxXxxxxXXXXXX* 

* INCREMENT MEAN ANOMALY & APPLY SECULAR VARIATION 

<^x*****rt*******x***rt**x*xx*x*x***xx*x****rt**xxx**xxxxxxx********* 

mnmot ion=xnow (1 ) * (-3/2) ; 

xnow ( 6) =twopi (xnow ( 6) +mnmot ion*Tloop) ; 

EA (k+1 ) =enewton (xnow (6) , xnow (6) , xnow (2) ) / 
cosTA= (cos (EA (k+1 ) ) -xnow (2 ) ) / (1 -xnow (2) *cos (EA (k-1 ) ) ) ; 
sinTA=sqrt ( 1 -xnow (2) *2 ) *sin (EA (k + 1 ) ) / (1- xnow (2) * cos (EA (k+1 ) ) ) ; 
r=xnow (1) * (1-xnow (2) *2) / (1+xnow (2) *cosTA) ; 

TA (k+1 ) =twopi (asin ( sinTA) ) ; 
if TA (k+1 ) <pi 
if cosTA<0 

TA (k+1 ) =pi-TA (k + 1 ) ; % Perform quadrant check 

end % for true anomaly, 

else 

if cosTA<0 

TA (k + 1 ) = 3 * p i - T A ( k i ) ; 
end 
end 

^-XXXXX*XX*XXXXXXXXXXXXXXXXXXXX'*XX*XXXXXXXXXXXXXXX')*-*'X?»'-**:XXXXXXXX'*'* 

* APPLY SECULAR VARIATIONS TO omega & OMEGA 

^★*x*xxxxx*xxxxxxxxxxXxxxx*xxx->rxxT*Xxxxx*xxxxxxxxxxxxxxx-*:xx*xxx*x* 

delTAsec=twopi (TA (k+1) -TA (k) ) ; 
semilat=xnow ( 1 ) * ( 1-xnow (2 ) *2 ) ; 

Q 1 = - 3 * J2 / s emi 1 a t * 2 ; 

SOMEGA=Ql * cos (xnow (3) ) / 2; 

Somega= (-Q1) * (2-5/2*sin (xnow (3) ) *2) / 2; 
de 1 sOMEGA=S OMEGA* delTAsec ; 
delsomega=Somega*delTAsec ; 
xnow (4 ) =xnow (4 ) +delsOMEGA; 
xnow (5) =xnow (5) +dels omega; 

^xx*****xx***x***x**x**xx***'**********x*xx*xx*'*xxxxx*xxxxx*x***'**k 

* CALCULATE ATMOSPHERIC DENSITY * 

f^*'**-****rtrt***rt******rtrtx*x*x**x'********x'**** , *x*x*xxxxx*xxxxx'***'*'!fc-* 

k2=3e- 6; 

k3=2 1/6378.135; 

atmdens=k2*exp ( (1-r) /k3) ; 

^************rt*X********rt*****rtrtrt******rt****rt*X***xrt************* 

* APPLY DRAG PERTURBATIONS TO a & e * 

<^******rt***rtrt*****rt***rt*******rt*******rt**rt*********rt**x********** 

kl = -4 * 637 8 . 135; 

deldrag= [kl*xnow (1) *2* (1 + xnow (2) ) *atmdens*delTAsec; 

kl*xnow (1) * (1-xnow (2) *2) *atmdens*delTAsec; 

0 ; 0 ; 0 ; 0 ] ; 
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xnow=xnow+ del drag; 

<^XXXXXX**X*************XXX**rt**’**rt*X******X*XXX*X*XXXXXXXXXXXXXXX 

* CALCULATE TRANSFORMATION MATRIX (ORBITAL TO INERTIAL) 

<£xxx********rt****rt******'x*xxxrtx*xx*xx**xxxxx*xxxxxxx*xxxxxxxxx*** 

cosperi = cos (xnow (5) ) ; 
sinperi=sin (xnow (5) ) ; 

cosascnd=cos (xnow ( 4 ) ) ; % Calculate sines & cosines of 

sinascnd=sin (xnow (4 ) ) ; % angular orbital elements. 

costilt=cos (xnow (3 ) ) ; 
sintilt=sin (xnow (3) ) ; 

Pl=cosperi*cosascnd-sinperi*sinascnd* costilt ; 

F2=-s inperi *cosascnd-cosperi *sinascnd*costilt; 
P4=cosperi*sinascnd+sinperi*cosascna*cost ilt ; 
P5=-sinperi*sinascna+cosperi x cosascnd*costilt; 
P7=sinperi*sintilt ; 

P£=cosperi * sintilt; 

S^**X*X****rt**X'*XX*XX*'*'**XX*'X* , **'** , 3k'rt'* , rt**'*XX'*XXTX*'XXXX>X>T*-XX****XX'* 

* UPDATE TIME * 

<^xxX*xx*******xx*****x*x**x*x*xx***x**x*xx**x*xxxxxxx**xxx**x*rt** 

Tgnwch=Tgnwch+Tloop; 

Tf ilter=Tf ilter+Tloop; 

5^XX7*'X*X**Xrt*'**x , *XX'*'X'*'*XX'*XXX'* - ****X*X'*'* , >k*X**'*xxx-*xxxxxxxXx-)*-'*xX**->*'?k- 

* CALCULATE OUT-OF-FENCE-PLANE COORDINATE OF SATELLITE POSITION * 

<^X*X***X****X'*Xrt , **'*'***XX'****XX*X'**'*XX*'*XX*XXXXXXXXXXXXXXXXXXXXXXX 

theta (k + 1 ) =twopi (earthrot *Tgnwch+lonx) ; 
sintheta=sin (theta (k+1) ) ; 
costheta-cos (theta (k+1) ) ; 
z f old=z fence ; 

zfence=.0031+(-sinfinc*costheta*Pl. . . 

-sinf inc * sintheta *P4 . . . 

+cosf inc*P7) *r*cosTA. . . % Calculate current magnitude 

+ (-sinfinc*costheta*P2 . . . % of out-of-plane coordinate, 
-sinf inc*sintheta*P5 . . . 

+cosfinc*P8) *r*sinTA; 

o 

* CHECK FOR FENCE PLANE INTERSECTION * 

O 

* NEGATIVE TO POSITIVE ?? * 

if zfence>0 
if zfoldcO 
n=n+l ; 
negtopos=n; 
xnow=xnow 1 ; 

[Tloop,Tfilter, Tgnwch, flag, xnow, r, cosTA, . . . 
sinTA, TA (k+1 ) , theta (k+1) , EA (k+1 ) , zfence] = . . . 
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inttemp (xnow, r, cosTA, . . . 

sinTA, TA (k + 1 ) , theca (k + 1 ) , EA (k + 1 ) , PI , P2 , ?4 , ... 

P5, P7 , P8 , Tgnwch, Tf ilter, Tloop) ; 
xnow=xnow ' ; 
end 
end 

* POSITIVE TO NEGATIVE ?? * 

£-*************************'***>'******x*-*******'*'****x************** 

O 

if zfence<0 
if zfold>0 
n=n + l ; 
postoneg=n; 
xnow=xnow 1 ; 

[Tloop, Tfilter, Tgnwch, flag, xnow, r, cosTA, . . . 
sinTA, TA(k+l) , theta (k+1) ,EA(k+l) , zfence]=. . . 
inttemp (xnow, r, cosTA, . . . 

sinTA, TA( k+1) , theta (k+1) , EA (k+1) , PI, P2, ?4, ... 

P5, P7, P8, Tgnwch, Tfilter, Tloop) ; 
xnow=xnow 1 ; 
end 
end 

TAold=TA (k+1 ) ; 

* SATELLITE IN NAVSPASUR WINDOW ?? - 

if flag==l 

* CALCULATE TIME FROM LAST FILTER UPDATE (LAST OBS) 

9'X**X********X**-*****X***r*****'****-***'*******-**'***-XXX*'*:*:*X******** 

C 



del obs=t imeobs-Tfi Iter ; 
while abs (delobs ) <1 
if length (observ (:, 1 )) >=p 
aelobs=t imeobs -Tfilter ; 

delT=delobs ; 

delTint (p) =delT*secperTU; 

Tint (p) = (timeobs*secperTU) /3600; 

Tupdate=t imeobs -Told; 
if TupdatecO 

Tupdate=Tupdate+TUperday; 

end 

Told=t imeobs ; 
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Tfilter=Tfilter+delT; 

Tgnwch=Tgnwch+delT; 



T X X * 



^X'*******'X***-**********'**'***'*** , ***'***********'*XX*XXXXXXXXXXXX 

* INCORPORATE PERIODIC VARIATIONS * 

^**X:X**********************-*****-***************'XX*xx**xxxxxx***t*-* 

xmean=xnow; 

daper=3 * J2*xnow ( 1 ) / (2 * semi 1 at A 2 * ( 1-xnow (2 ) A 2 ) ) * . . . 

sin (xnow (3) ) A 2*cos (2* (xnow (5) +TA (k+1 ) ) ) ; 
deper=3 * J2 / (2*semilat A 2) * ( (l-3/2*sin(xnow(3) ) A 2 ) * cosTA+ . . . 
l/4*sin(xnow(3) ) A 2 *cos (2* xnow (5) +TA (k+1 ) ) + . . . 

7/1 2* sin (xnow (3) ) A 2*cos (2* xnow (5) +3*TA (k+1) ) ) ; 
diper=-3* J2/ ( 8 * semi 1 at A 2 ) *sin (2*xnow (3) ) * . . . 

cos (2* (xnow (5) +TA (k+1 ) ) ) ; 
dOMper=3*J2/ (4 * semi lat A 2 ) *cos (xnow (3) ) * . . . 
sin (2* (xnow (5) +TA(k+l) ) ) ; 

dnewper=-3* J2/ ( 4 * semi 1 at A 2 ) * ( 1 -5/2 * sin (xnow ( 3 ) ) A 2 ) . . . 

’'sin (2* (xnow (5) +TA(k+l) ) ) ; 
dtemp=3*J2/ (2 * semi 1 at A 2 ) * (1 /xnow (2) * . . . 

(l-3/2*sin (xnow (3) ) A 2) *sinTA- . . . 

1 /4 * sin (xnow (3) ) A 2*sin(2* xnow (5)-rTA(k + i))-r... 

7/12*sin (xnow (3) ) A 2 * sin (2* xnow (5) + 3*TA(k-rl) ) ) -r . . . 

1/2* ( (1-3/2* sin (xnow (3) ) A 2 ) * sin (2 * TA (k+ 1 ) ) - . . . 

( 1 -5/2 * sin (xnow ( 3 ) ) A 2 ) * sin (2 * (xnow (5 ) +TA (k + 1 ) ) ) + . . . 

3/8 * sin (xnow ( 3 ) ) A 2*sin(2*xnow(5)+4*TA(k + l) ) ) ) ; 
domper=dtemp; 
dMAper=dnewper-domper ; 

varper= [daper ; deper; aiper ; dOMper ; domper ; dMAper ] ; 
xnow=xnow+varper; xnow (6) =twopi (xnow (6) ) ; 

EA (k+1 ) =enewton (xnow ( 6) , xnow ( 6) , xnow (2 ) ) ; 
cosTA= (cos (EA(k+l) ) -xnow (2) ) / (1-xnow (2) *cos (EA (k+1) ) ) ; 
sinTA=sqrt (l-xnow(2) A 2) * sin (EA (k-t-1 ) ) / (l-xnow(2) ’cos (EA(k + l) ) ) / 
r=xnow ( 1 ) * ( 1-xnow (2 ) A 2 ) / ( 1 + xnow ( 2 ) * cosTA) ; 

TA (k+1 ) =twopi (as in (sinTA) ) ; 
if TA (k+1 ) <pi 
if cosTA<0 
TA (k + 1 ) =pi-TA (k + 1 ) ; 
end 
else 

if cosTA<0 

TA (k + 1 ) =3*pi-TA (k + 1 ) ; 
end 
end 



% Perform quadrant check 
% for true anomaly. 



* CALCULATE ANGULAR DISPLACEMENT OF RECEIVER SITE FROM I -AXIS * 
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phi=twopi (earthrot *Tgnwch+lon (observ (p, 3) ) ) ; 

sinphi=sin (phi) ; 

sinlat=sin (lac (observ (p, 3) ) ) ; 
coslat = cos (lat (observ (p, 3) ) ) ; 
sinlon=sin (Ion (observ (p, 3) ) ) ; 
coslon=cos (Ion (observ (p, 3) ) ) / 
cosperi=cos (xnow ( 5) ) ; 
sinperi = sin (xnow (5) ) / 
cosascnd=cos (xnow (4 ) ) ; 
sinascna=sin (xnow (4) ) ; 
cost i lt=cos (xnow (3) ) ; 
sintilt=sin (xnow (3) ) ; 

Pi=cosperi *cosascnd-sinperi*sinascnd*costilt ; 

?2 = -sinperi*cosascnd-cosperi , 'sinascnd’ : costilt ; 

F4=ccsperi , ■ s inascnd+sinper i * cosascnd* costilt ; 

P5 = -sinperi ’'sinascnd + cosperi* cosas end* costilt; 



P7 = sinperi* sint ilt ; 



P8=cosperi*sintilt; 

sinaz=sin (az (observ (p, 3) ) ) / 
cosaz = cos (az (observ (p, 3) ) ) ; 

Hi =c os lat *cosphi ; 

H2=coslat *sinphi; 

H3=sinlat ; 

H4 = -sinphi *cosaz-sinlat * cosphi *sinaz ; 
H5 = cosaz*cosphi-sinaz * sin lat *sinphi ; 
H6=sinaz*coslat; 

H7 = sinaz*sinphi-cosaz* sin lat * cosphi / 
H8 = -sinaz * cosphi -cos az*sinlat *sinphi; 
H9=cosaz*coslat ; 

R1=H4 *P1 + H5*P4 + H6*P7; 
R2=H4*P2+H5*P5+H6*P8; 

R4=H7 *P1 +H8 *P4+H9*P7 / 
R5=H7*P2+H8*P5+H9*P8; 
R7=H1*P1+H2*P4+H3*P7; 

R8=H1 *P2 + H2 *P5+H3*P8; 



o 
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* ESTIMATE MEASUREMENT * 

S£*X*-******-************-****-***X-**********'*'***'******'**-**-*> r^r^Tr****** 

Rsite=sqrt (1- (2-1/298.26) /298.26*sinlat A 2) +alt (cfcserv(p, 3) ) ; 
rof f set = . 0031 ; 

rhoH=r* (R7*cosTA+R8*sinTA) -Rsite; 

rhoE=r* (Rl*cosTA+R2*sinTA) ; 

rhoN=r* (R4 *cosTA+R5*sinTA) +roffset; 

rho=sqrt (rhoH A 2+rhoE A 2+rhoN A 2) ; 

zest ( : , p) = [ rhoE/ rho; rhoN/rho] ; 

size=xnow (1) ; shape=xnow (2) ;tilt=xnow (3) ; 
ascnd=xnow ( 4 ) ; per i=xnow (5 ) ; MA=xnow ( 6) ; 

* CALCULATE COVARIANCE MATRIX P * 

O 

J2=l . 082 64 5e -3; 
sinEA=sin (EA (k+1 ) ) ; 
cosEA=cos (EA (k+1 ) ) ; 

sindecl=s intilt* sin (TA (k+1 ) +peri) ; 
semilat=size* (l-shape A 2) ; 

G 

* CALCULATE APPROXIMATE delta (TRUE ANOMALY ) 

G 

mnmot ion=xnow ( 1 ) A ( - 3/2 ) ; 

delMA=mnmotion*Tupdate; 

delTA=delMA; 

Alphi=- (l+shape A 2) *sinEA+2*shape*sinEA*cosEA; 

A2phi=l- shape* cosEA; 

A3phi = 2-5 /2 * sintilt A 2; 

A4phi=l-3*sindecl A 2; 

A5phi=l -shape A 2; 

A6phi=sqrt (l/size+J2/r A 3*A4phi) ; 

A7phi=l+shape* cosTA; 

dOMEGAda=3 * J2 / (size*semilat A 2) *costilt *delTA; 
dOMEGAde=- 6* J2 * size* shape/ semi lat A 3*cos tilt *delTA; 
dOMEGAdi=3* J2/ (2*semilat A 2 ) *sintilt*delTA; 

dOMEGAdM=-3* J2/ (2*semilat A 2 ) *costilt *Alphi/ (sinTA*A2phi A 3) ; 

domegada=-3* J2/ (size*semilat A 2) *A3phi*delTA; 
domegade=6* J2 *size*shape/semilat A 3*A3phi*delTA; 
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domegadi = -15* J2/ (2*semilat~2) * sintilt * cost ilt * del TA; 
domegadM=-3* J2/ (2*semilat~2) *A3phi*Alphi/ ( sinTA* A2phi A 3 ) ; 

drkdak=A5phi /A7phi ; 

dMAda=3/2*Tupdate*A6phi* (-l/size~2+J2*A4phi* (-3/r A 4) *crkdak) ; 

drkdek= (-2*size*shape* (l-shape*cosTA) -semilat* cosTA) /A7phi~2; 
dMAde=3/2*Tupdate*A6phi* J2*A4phi* (-3/r~4) *drkdek; 

dsindkdi=sin (TA (k+1) +peri) *costilt; 
dMAdi=3*Tupdate*A6phi* (-3*J2/r^3) *sindecl*dsindkdi; 

asindkdo=sintilt*cos (TA (k+1) +peri) / 
dMAdomega=3*Tupdate*A6phi* (-3*J2/r~3) *sindecl*dsindkdo; 
drdTA=semilat*shape*sinTA/A7phi / '2; 
dsinddTA=dsindkdo; 

dmessdTA=-3* J2/r^2*drdTA- (-9* J2/r^4*sindecl‘'2*drdTA+ . . . 

6* J2/r^3*sindecl*dsinddTA) ; 
dTAdMA=-Alphi / (sinTA*A2phi A 3) / 
dMAdM=l+3/2*Tupdate*A6phi*dmessdTA*dTAdMA; 

dada=l + 2*kl ’'size* (1 + shape) *atmaens*delTA; 

dade=kl*size~2*atmdens*delTA; 

dadM=-kl *size^2* (1+shape) * atmdens*dTAdMA; 

deaa=kl *A5phi *atmdens*delTA; 
dede-l-2*kl *size* shape* atmdens *delTA; 
dedM=-kl * size*A5phi * atmdens *dTAdMA; 

phimat=[dada dade 0 0 0 dad.M; 
deda dede 0 0 0 dedN; 

001000; 

dOMEGAda dOMEGAde dOMEGAdi 1 0 dOMEGAdK; 
domegada domegade domegaai 0 1 domegadM; 
dMAda dMAde dMAdi 0 dMAdomega dMAdM] ; 

c 

* CALCULATE PLANT NOISE * 

J2 = l . 0 82 6e- 3 ; 

semilat=xnow (1) * (1-xnow (2) *2) ; 

siga2= (3* J2*xnow (1) / (2*semilat~2* (1-xnow (2) ~2) ) * . . . 

sin (xnow (3) ) ^2*003 (2* (xnow (5) +TA (k+1) ) ) ) ~2; 
sige2=(3*J2/ (2 * semi lat ~2 ) * ( (l-3/2*sin (xnow (3) ) ~2) *cosTA+ . . . 
sigi2= (3*J2/ (8* semilat ~2 ) *sin (2*xnow(3) ) * . . . 

cos (2* (xnow (5) +TA (k+1) ) ) ) A 2; 
mnmotion=xnow ( 1 ) A (-3/2 ) ; 

MApert=mnmot ion* 3* J2/ (2* semilat ^2) * (l-3/2*sin(xnow(3) ) ~2) * . . . 
sqrt (1-xnow (2)^2) ; 

sigOMEG2= (3* J2/ (2* semilat *2) *mnmot ion* cos (xnow (3) ) *Tupdate) ~2; 
sigomeg2= (3*J2/ (2*semilat~2) *mnmotion* . . . 



88 



(2-5/2 * sin (xnow (3) ) "2) *Tupdate) ~2; 
sigMA2= (MApert *Tupdate) ~2; 

Q= [ le-2 * siga2 0 0 0 0 0 
0 le-3*sige2 0000 
0 0 sigi2 000 
000 sigOMEG2 0 0 
0000 sigomeg2 0 
00000 sigMA2 ] . *le-3; 

S-X*'*X'** , *X'*'*'*'*'*X'*'***'*'*X**X'**')kx*x*'**X*'*XxXXXXXXXXXx*-*-j»-**'*'*'*'X*XXxX , *X 

o 

* PROPAGATE ERROR COVARIANCE * 

Pkk=Pkkml ; 

Pkkml=phimat *Pkk*phimat ’ +Q; 

S_**X*X*XX****XX****xX*XXX***'*X*X*XXXXXXxxXX**XX'**'yX:X:'******r**x***X 

o 

- CALL IN OBSERVATION DATA * 

S-XxxX- , XXXXX*X'*XXXXxX')k*X*XXXxXXXXXXXXXXXXXX*XXxXXx-*Xx*xxxXXXX'XXXXXX 

c 

cosalpha=sin (observ (p, 4) *pi/180) ; 
cosbeta=sin (observ (p, 5) *pi/180) / 
z ( : ,p) = [cosalpha; cosbeta] ; 

S^XxxxXXXX*XXXXXX*****'*XXXXxXX'*'*X*XX*X*>k-'*XXXxXxxXxXXxXXXX*:XX*:XX*')k* 

* CALCULATE LINEARIZED MEASUREMENT MATRIX H * 

S^xx*XX***X**X*XXxXXXXXX*XXXXXx*X*X*xXxXxxXXXXXXXX'*'XXX*7'rxxxXXx***X 

c 

a 1=1 -xnow (2) ^2; 

a 2=1+ xnow (2) * cosTA; 

a3=2*xnow(2) + cosTA* ( 1 + xnow (2 ) * 2 ) ; 

a4=2*xnow (2) *sinEA*cosEA; 

a5= 1+xnow (2 ) ^2 / 

a6=-a5* sinEA+a4 ; 

a7=l -xnow (2 ) * cosEA; 

drda=al /a2; 

drde=-xnow (1) *a3/a2*2; 
dcosTAdT=-sinTA; 

dTAdEA= (-a5*sinEA+a4 ) / (-sinTA*a7~2 ) ; 
dE AdMA= 1 / a 7 ; 

dc O S T AdM= d COS T AdT * dT AdE A * dE AdMA ; 
dsinTAdT= cosTA; 

d s i n TAdM=ds inTAdT * dTAdEA * dEAdMA / 

drdTA=xnow ( 1 ) *xnow (2 ) *al *sinTA/a2 , '2 ; 
drdMA=drdTA* dTAdEA* dEAdMA; 

dRldi=sinperi*sinascnd*sintilt *H4- . . . 

sinperi*cosascnd*sintilt*H5+ . . . 
sinperi* costilt *H6; 
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dR2di=cosperi * sinascnd* sintilt *H4 - . . . 

cosperi * cosascnd* sint ilt *K5-r . . . 
cosperi* costilt *H6; 
dR4di=sinperi*sinascnd*sintilt *K7- . . . 

s inperi* cosascnd* sint ilt *H8+ . 
sinperi*costilt*H9; 
dR5di = cosperi* sinascnd* sintilt *K7- . 

cosperi* cosascnd* sintilt *H6+ . 
cosperi * cost ilt *H9; 
dR7di=s inperi* sinascnd* sint ilt *H1 - . 

s inperi* cosascnd* sint i It *H2+ . 
s inperi * cost ilt *H3; 
dR8di = cosperi * sinascnd* sintilt *H1- . 

cosperi* cosascnd* sint i It *K2+ . 
cosperi * cost ilt *H3; 

dRldOM=-P4*H4+Pl*H5; 
dR2dOM=-P5*H4+P2*H5; 
dR4aOM=-P4*H7+Pl*H8; 
dR5dOM=-P5*H7+P2*H8; 
dR7 dOK= -P 4 * HI +P 1 * H2 ; 
dR8dOR=-?5*Hl+P2*K2; 

dRldom=P2*H4+P5*H5+P8*H6; 
dR2dom=-Pl*H4-P4*K5-P7*H6; 
dR4aom=P2 *H7+P5*H6+P8*HS; 
dR5dom=-Pl*H7-P4*H8-P7*H9; 
dR7dom=P2*Kl+P5*H2+P8*H3; 
dR8dom=-P 1 *H1-P4*H2-P7*H3; 



dpHda= (R7 * cosTA+RS * sinTA) *drda; 
dpEda= (Rl*cosTA+R2*sinTA) *drda; 
dpNda= (R4 *cosTA+R5*sinTA) *drda; 
dpda= (rhoH*dpHda+rhoE*dpEda+rhoN*dpNda) /rho 

dpHde= (R7*cosTA+R8*sinTA) *drde; 
dpEde= (R1 *cosTA+R2 * sinTA) *drae; 
dpNde= (R4 *cosTA+R5*sinTA) *drde; 
dpde= (rhoH*dpHde+rhoE*dpEde+rhoN*dpNde) /rho 

dpHdi= (dR7di *cosTA+dR8di * sinTA) *r; 
dpEdi= (aRldi*cosTA+dR2ai*sinTA) *r; 
dpNdi= (aR4di*cosTA+dR5di*sinTA) *r; 
dpdi= (rhoH*dpHdi+rhoE*dpEdi+rhoN*dpNdi) /rho 

dpHdOM= ( dR7 dOM *c os TA+dR8dOM* sinTA) *r; 
dpEdOM= (dRldOM*cosTA+dR2dOM* sinTA) *r; 
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dpNdOM= (dR4dOM*ccsTA+dR5dOM*sinTA) *r; 

dpdOM= (rhoH*dpHdOM+rhoE*dpEdOM+rhoN*dpNdOM) /rho; 

dpHdom= (dR7dom* cosTA+dR8dom* sinTA) * r ; 
dpEdom= (dRldom*cosTA+dR2dom* sinTA) *r; 
dpNdom= (dR4dom*cosTA+dR5dom* sinTA) *r; 
dpdom= (rhoH*dpHdom+rhoE*dpEdom+rhoN*dpNdom) /rho; 

dpHdMA= (R7*dcosTAdM+R8*dsinTAdM) *r+ (R7 *cosTA+R8 * 
sinTA) *drdMA; 

dpEdMA= (R1 *dcosTAdM+R2 *dsinTAdM) *r+ (R1 *cosTA-rR2 * 
sinTA) *drdMA; 

dpNdMA= (R4 *dcosTAdM+R5*dsinTAdM) *r+ (R4 *ccsTA+R5* 
sinTA) *drdMA; 

dpdKA= (rhoH*dpHdMA+rhoE*dpEdMA+rhoN*dpNdMA) /rho; 

drhoEda=dpEda/rho-rhoE*dpda/rho~2 ; 
drhoNda=dpNda/ rho-rhoN*dpda/rho A 2 ; 

drhoEde=apEde/rho-rhoE*dpde/rhc~2; 
drhoNde=dpNde/rho-rhoN*dpde/rho~2 ; 

drhoEdi=dpEdi/rho-rhoE*dpdi/rho~2; 
drhoNdi=dpNdi /rho-rhoN , 'dpdi/rho~2 ; 

drhoEdOM=dpEdOM/rho-rhoE*dpdOM/rho^2; 
drhoNdOM=dpNdOM/rho-rhoN*dpdOM/rho A 2 ; 

drhoEdom=dpEdom/rho-rhoE*dpdom/rho~2 ; 
drhoNdom=dpNdom/rho-rhoN*dpdom/rho A 2 ; 

drhoEdMA=dpEdMA/rho-rhoE*dpdMA/rho^2 ; 
arhoNdMA=dpNdMA/rho-rhoN* dpdMA/rho A 2 ; 

H=[drhoEda drhoEde arhoEdi drhoEaOM arhoEdom drhoEdMA; 
drhoNda drhoNde drhoNdi drhoNdOM drhoNdom drhoNdMA] ; 

* CALCULATE KALMAN GAIN * 

G=Pkkml*H' *inv (H*Pkkml *H ' +R) ; 
res id ( : , p) =z ( : , p) -zest ( : , p) ; 



%**************************************************************** 
* CALCULATE IMPROVED COVARIANCE * 
%**************************************************************** 
Pkk= (eye (6) -G*H) *Pkkml; 
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Pkeep(:,p*6+1:6* (p+1 ) ) =Pkk; 

* CALCULATE IMPROVED MEAN ORBITAL ELEMENTS * 

^*X**XX**XX*X**X , *X*XX*X**XXXXX**XX*X'*XX*XXXXX*X:I*XXXXXXX:*'*»:**X*X** 

xnow=xmean ; 

G* res id ( : , p) ; 
xnow=xnow+G*resid ( : , p) ; 
xnow ( 6) =twopi (xnow ( 6) ) ; 



xint ( 1 : 2 , p) =xnow (1:2) ; 

xint (3 : 6, p) =xnow (3:6) *18 0/pi; 

XX*X**X***X*-*X , '******X**X****X*XX*-*-**XXXXXXX*XXXx:*XXXXXXXxXXXXXX* 

* CALCULATE IMPROVED ECCENTRIC ANOMALY * 

■*-XXXX*'****XX*XX*X**XXXXXX*XX*XX*XX*XXX**X*XXXXXXXX'*XXXXX*XXXXXXX 

EA (k + 1 ) =enewton (xnow ( 6) , xnow ( 6) , xnow ( 2 ) ) ; 

■F>--^7irXXX^X>t>kr*^-^')trxXX'*^**X**T^*'*' , ^y^XXXXX^^XXXXX'*'3*:XXX'5»X'?tXXXX'J»'*^XXX>r^: , >: 

* SOLVE FOR IMPROVED cos (TA) , sin(TA), r, sin(decl), n * 

■F>-*X'FXT(tXX'^XX7»X-^-^^'XXX*-*-X**XXXXXX^*XXXXXXXX'J*-l*-X , J»XX7»X7*>t-j»'F7tTrxXXXXX'*'X>k- 

cosTA= (cos (EA(k + l) ) -xnow (2) ) / (1-x.now (2) *cos (EA (k-rl) ) ) ; 
sinTA=sqrt (1- xnow (2) ^2) *sin (EA (k+ 1 ) ) / ( i-xnow (2 ) . . . 

*cos (EA (k+1) ) ) ; 

r=xnow ( 1 ) ’ ( 1-xnow (2 ) *'2 ) / ( 1 -rxnow (2 ) * cosTA; ; 



TA (k-1 ) =twopi (as in ( sinTA) ) ; 
if TA (k+1 ) <pi 
if cosTA<0 

TA (k+1 ) =pi-TA (k+1 ) ; % Perform quadrant check 

end % for true anomaly, 

else 

if cosTA<0 

TA (k + 1 ) =3*pi-TA (k + 1 ) ; 
end 
end 



TAold=TA (k+1) ; 









XXxXXX-*X*XXXXXXXX*XX 



if length (observ (:, 1 )) >p 
p=p+l; 

time=observ (p, 2) ; 
hrs=fix (time*le-4) ; 
min=fix (rem (time, le4) * le-2) ; 
sec=rem (time, le2) ; 

timeobs= (hrs*3600+min*60+sec) /secperTU; 
delobs=t imeobs-Tf ilter; 
else 

aelobs=2 ; 
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